<!-- saved from url=(0014)about:internet -->
<pre class="code">
<span class="srcline"><span class="lineno"><a href="1,1" id="srcline1">  1</a></span><span class="line"><span class="keyword">function</span> [<span class="var type1" id="S2T1U3">tout</span>, <span class="var type1" id="S3T35U4">xout</span>, <span class="var type1" id="S4T1U5">ErrorFlag</span>] = ode78_(<span class="var type1" id="S5T3U8">tspan</span> , <span class="var type1" id="S6T9U9">x0</span> , <span class="var type1" id="S7T4U10">auxdata</span>)</span></span>
<span class="srcline"><span class="lineno"><a href="1,2" id="srcline2">  2</a></span><span class="line"><span class="comment">%#codegen</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,3" id="srcline3">  3</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,4" id="srcline4">  4</a></span><span class="line"><span class="comment">% ------------------------------------------------------------------</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,5" id="srcline5">  5</a></span><span class="line"><span class="comment">% rk78 integration</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,6" id="srcline6">  6</a></span><span class="line"><span class="comment">%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,7" id="srcline7">  7</a></span><span class="line"><span class="comment">% [tout, xout, ErrorFlag] = ode78_chen_(tspan , x0 , auxdata)</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,8" id="srcline8">  8</a></span><span class="line"><span class="comment">%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,9" id="srcline9">  9</a></span><span class="line"><span class="comment">% Input argumuents:</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,10" id="srcline10"> 10</a></span><span class="line"><span class="comment">% -------------------------------------------------------------------------</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,11" id="srcline11"> 11</a></span><span class="line"><span class="comment">% tspan            [1x2 double]           time array                                  [ TU ]        </span></span></span>
<span class="srcline"><span class="lineno"><a href="1,12" id="srcline12"> 12</a></span><span class="line"><span class="comment">% x0                 [14x1 double]         state and costate variables     [  -  ]</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,13" id="srcline13"> 13</a></span><span class="line"><span class="comment">%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,14" id="srcline14"> 14</a></span><span class="line"><span class="comment">% Output argumuents:</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,15" id="srcline15"> 15</a></span><span class="line"><span class="comment">% -------------------------------------------------------------------------</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,16" id="srcline16"> 16</a></span><span class="line"><span class="comment">% tout              [1x1 double]            position vector          [  -  ]</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,17" id="srcline17"> 17</a></span><span class="line"><span class="comment">% xout             [14x1 double]          velocity vector           [  -  ]</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,18" id="srcline18"> 18</a></span><span class="line"><span class="comment">% ErrorFlag   [1x1 double]             error flag                   [  -  ]</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,19" id="srcline19"> 19</a></span><span class="line"><span class="comment">%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,20" id="srcline20"> 20</a></span><span class="line"><span class="comment">% External functions called:</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,21" id="srcline21"> 21</a></span><span class="line"><span class="comment">% -------------------------------------------------------------------------</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,22" id="srcline22"> 22</a></span><span class="line"><span class="comment">% fy_.m</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,23" id="srcline23"> 23</a></span><span class="line"><span class="comment">%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,24" id="srcline24"> 24</a></span><span class="line"><span class="comment">% Copyright (C) 2015/06/22 by Chen Zhang, </span></span></span>
<span class="srcline"><span class="lineno"><a href="1,25" id="srcline25"> 25</a></span><span class="line"><span class="comment">% BeiHang University, chenzhang.buaa@gmail.com</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,26" id="srcline26"> 26</a></span><span class="line"><span class="comment">% ------------------------------------------------------------------</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,27" id="srcline27"> 27</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,28" id="srcline28"> 28</a></span><span class="line"><span class="comment">% The coefficients of rk78</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,29" id="srcline29"> 29</a></span><span class="line"><span class="mxinfo" id="T26:U7"><span class="var type1" id="S8T26U13">c_i</span>=  <span class="mxinfo" id="T26:U9">[ 1/18, 1/12, 1/8, 5/16, 3/8, 59/400, 93/200, 5490023248/9719169821, 13/20, 1201146811/1299019798, 1, 1]'</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,30" id="srcline30"> 30</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="var type1" id="S9T27U51">a_i_j</span> = <span class="mxinfo" id="T27:U12">[ 1/18, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,31" id="srcline31"> 31</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    1/48, 1/16, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,32" id="srcline32"> 32</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    1/32, 0, 3/32, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,33" id="srcline33"> 33</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    5/16, 0, -75/64, 75/64, 0, 0, 0, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,34" id="srcline34"> 34</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    3/80, 0, 0, 3/16, 3/20, 0, 0, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,35" id="srcline35"> 35</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    29443841/614563906, 0, 0, 77736538/692538347, -28693883/1125000000, 23124283/1800000000, 0, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,36" id="srcline36"> 36</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    16016141/946692911, 0, 0, 61564180/158732637, 22789713/633445777, 545815736/2771057229, -180193667/1043307555, 0, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,37" id="srcline37"> 37</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    39632708/573591083, 0, 0, -433636366/683701615, -421739975/2616292301, 100302831/723423059, 790204164/839813087, 800635310/3783071287, 0, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,38" id="srcline38"> 38</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    246121993/1340847787, 0, 0, -37695042795/15268766246, -309121744/1061227803, -12992083/490766935, 6005943493/2108947869, 393006217/1396673457, 123872331/1001029789, 0, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,39" id="srcline39"> 39</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    -1028468189/846180014, 0, 0, 8478235783/508512852, 1311729495/1432422823, -10304129995/1701304382, -48777925059/3047939560, 15336726248/1032824649, -45442868181/3398467696, 3065993473/597172653, 0, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,40" id="srcline40"> 40</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    185892177/718116043, 0, 0, -3185094517/667107341, -477755414/1098053517, -703635378/230739211, 5731566787/1027545527, 5232866602/850066563, -4093664535/808688257, 3962137247/1805957418, 65686358/487910083, 0, 0;</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,41" id="srcline41"> 41</a></span><span class="line"><span class="mxinfo" id="T27:U10"><span class="mxinfo" id="T27:U12">    403863854/491063109, 0, 0, -5068492393/434740067, -411421997/543043805, 652783627/914296604, 11173962825/925320556, -13158990841/6184727034, 3936647629/1978049680, -160528059/685178525, 248638103/1413531060, 0, 0]'</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,42" id="srcline42"> 42</a></span><span class="line"><span class="mxinfo" id="T28:U13"><span class="var type1" id="S10T28U362">b_8</span> = <span class="mxinfo" id="T28:U15">[ 14005451/335480064, 0, 0, 0, 0, -59238493/1068277825, 181606767/758867731,   561292985/797845732,   -1041891430/1371343529,  760417239/1151165299, 118820643/751138087, -528747749/2220607170,  1/4]'</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,43" id="srcline43"> 43</a></span><span class="line"><span class="mxinfo" id="T28:U16"><span class="var type1" id="S11T28U402">b_7</span> = <span class="mxinfo" id="T28:U18">[ 13451932/455176623, 0, 0, 0, 0, -808719846/976000145, 1757004468/5645159321, 656045339/265891186,   -3867574721/1518517206,   465885868/322736535,  53011238/667516719,                  2/45,    0]'</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,44" id="srcline44"> 44</a></span><span class="line"><span class="mxinfo" id="T1:U19"><span class="var type1" id="S12T1U439">pow</span> = <span class="mxinfo" id="T1:U21">1/8</span></span>; <span class="comment">% power for step control</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,45" id="srcline45"> 45</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,46" id="srcline46"> 46</a></span><span class="line"><span class="comment">% Maximum step size</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,47" id="srcline47"> 47</a></span><span class="line"><span class="mxinfo" id="T1:U22"><span class="var type1" id="S13T1U445">hmax</span> = <span class="mxinfo" id="T1:U24">(<span class="mxinfo" id="T1:U25"><span class="mxinfo" id="T1:U26"><span class="var type1" id="S5T3U450">tspan</span>(<span class="mxinfo" id="T1:U28">2</span>)</span> - <span class="mxinfo" id="T1:U29"><span class="var type1" id="S5T3U453">tspan</span>(<span class="mxinfo" id="T1:U31">1</span>)</span></span>)/<span class="mxinfo" id="T1:U32">2.5</span></span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,48" id="srcline48"> 48</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,49" id="srcline49"> 49</a></span><span class="line"><span class="comment">% Initial step size</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,50" id="srcline50"> 50</a></span><span class="line"><span class="mxinfo" id="T1:U33"><span class="var type1" id="S14T1U458">h</span> = <span class="var type1" id="S13T1U459">hmax</span></span> ; <span class="keyword">if</span> <span class="mxinfo" id="T2:U36"><span class="var type1" id="S14T1U463">h</span> &gt; <span class="mxinfo" id="T1:U38">0.5</span></span> ; <span class="mxinfo" id="T1:U39"><span class="var type1" id="S14T1U467">h</span> = <span class="mxinfo" id="T1:U41">0.5</span></span> ; <span class="keyword">end</span> ; <span class="keyword">if</span> <span class="mxinfo" id="T2:U42"><span class="var type1" id="S14T1U472">h</span> &gt; <span class="var type1" id="S13T1U473">hmax</span></span> ; <span class="mxinfo" id="T1:U45"><span class="var type1" id="S14T1U476">h</span> = <span class="var type1" id="S13T1U477">hmax</span></span> ; <span class="keyword">end</span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,51" id="srcline51"> 51</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,52" id="srcline52"> 52</a></span><span class="line"><span class="comment">% Relative error tolerance (applies to all components of the solution vector)</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,53" id="srcline53"> 53</a></span><span class="line"><span class="mxinfo" id="T1:U48"><span class="var type1" id="S15T1U480">tol</span> = <span class="mxinfo" id="T1:U50"><span class="var type1" id="S7T4U482">auxdata</span>.tol</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,54" id="srcline54"> 54</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,55" id="srcline55"> 55</a></span><span class="line"><span class="mxinfo" id="T1:U52"><span class="var type1" id="S16T1U486">t0</span> = <span class="mxinfo" id="T1:U54"><span class="var type1" id="S5T3U488">tspan</span>(<span class="mxinfo" id="T1:U56">1</span>)</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,56" id="srcline56"> 56</a></span><span class="line"><span class="mxinfo" id="T1:U57"><span class="var type1" id="S17T1U492">tfinal</span> = <span class="mxinfo" id="T1:U59"><span class="var type1" id="S5T3U494">tspan</span>(<span class="mxinfo" id="T1:U61">2</span>)</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,57" id="srcline57"> 57</a></span><span class="line"><span class="mxinfo" id="T1:U62"><span class="var type1" id="S18T1U498">t</span> = <span class="var type1" id="S16T1U499">t0</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,58" id="srcline58"> 58</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,59" id="srcline59"> 59</a></span><span class="line"><span class="comment">% Minimum step size</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,60" id="srcline60"> 60</a></span><span class="line"><span class="mxinfo" id="T1:U65"><span class="var type1" id="S19T1U502">hmin</span> = <span class="mxinfo" id="T1:U67"><span class="mxinfo" id="T1:U68">16*eps</span>*<span class="mxinfo" id="T1:U69">abs(<span class="var type1" id="S18T1U510">t</span>)</span></span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,61" id="srcline61"> 61</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,62" id="srcline62"> 62</a></span><span class="line"><span class="comment">% Initialization</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,63" id="srcline63"> 63</a></span><span class="line"><span class="mxinfo" id="T1:U71"><span class="var type1" id="S22T1U513">n_reject</span> = <span class="mxinfo" id="T1:U73">0</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,64" id="srcline64"> 64</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,65" id="srcline65"> 65</a></span><span class="line"><span class="comment">% constant for step rejection</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,66" id="srcline66"> 66</a></span><span class="line"><span class="mxinfo" id="T1:U74"><span class="var type1" id="S23T1U517">reject</span> = <span class="mxinfo" id="T1:U76">0</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,67" id="srcline67"> 67</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,68" id="srcline68"> 68</a></span><span class="line"><span class="comment">% counter for function evaluations</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,69" id="srcline69"> 69</a></span><span class="line"><span class="mxinfo" id="T1:U77"><span class="var type1" id="S24T1U521">fk</span> = <span class="mxinfo" id="T1:U79">0</span></span> ; <span class="mxinfo" id="T1:U80"><span class="var type1" id="S4T1U525">ErrorFlag</span> = <span class="mxinfo" id="T1:U82">0</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,70" id="srcline70"> 70</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,71" id="srcline71"> 71</a></span><span class="line"><span class="mxinfo" id="T9:U83"><span class="var type1" id="S25T9U529">x</span> = <span class="mxinfo" id="T9:U85"><span class="var type1" id="S6T9U531">x0</span>(:)</span></span> ; <span class="comment">% start point</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,72" id="srcline72"> 72</a></span><span class="line"><span class="mxinfo" id="T11:U87"><span class="var type1" id="S26T11U535">f</span> = <span class="mxinfo" id="T11:U89"><span class="var type1" id="S25T9U537">x</span>*<span class="mxinfo" id="T10:U91">zeros(1,13)</span></span></span> ;  <span class="comment">% array f for RHS calculation</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,73" id="srcline73"> 73</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,74" id="srcline74"> 74</a></span><span class="line"><span class="mxinfo" id="T1:U92"><span class="var type1" id="S2T1U544">tout</span> = <span class="var type1" id="S18T1U545">t</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,75" id="srcline75"> 75</a></span><span class="line"><span class="mxinfo" id="T35:U95"><span class="var type1" id="S3T35U548">xout</span> = <span class="mxinfo" id="T35:U97"><span class="var type1" id="S25T9U550">x</span>.'</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,76" id="srcline76"> 76</a></span><span class="line"><span class="comment">% tau = tol * max(norm(x,'inf'), 1) ;  % accuracy</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,77" id="srcline77"> 77</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,78" id="srcline78"> 78</a></span><span class="line"><span class="comment">% The main loop</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,79" id="srcline79"> 79</a></span><span class="line"><span class="comment">% count = 1;</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,80" id="srcline80"> 80</a></span><span class="line"><span class="keyword">while</span> (<span class="mxinfo" id="T2:U99"><span class="var type1" id="S18T1U555">t</span> &lt; <span class="var type1" id="S17T1U556">tfinal</span></span>) &amp;&amp; (<span class="mxinfo" id="T2:U102"><span class="var type1" id="S14T1U559">h</span> &gt;= <span class="var type1" id="S19T1U560">hmin</span></span>)</span></span>
<span class="srcline"><span class="lineno"><a href="1,81" id="srcline81"> 81</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,82" id="srcline82"> 82</a></span><span class="line">    <span class="keyword">if</span> <span class="mxinfo" id="T2:U105">(<span class="mxinfo" id="T1:U106"><span class="var type1" id="S18T1U566">t</span> + <span class="var type1" id="S14T1U567">h</span></span>) &gt; <span class="var type1" id="S17T1U568">tfinal</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,83" id="srcline83"> 83</a></span><span class="line">        <span class="mxinfo" id="T1:U110"><span class="var type1" id="S14T1U571">h</span> = <span class="mxinfo" id="T1:U112"><span class="var type1" id="S17T1U573">tfinal</span> - <span class="var type1" id="S18T1U574">t</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,84" id="srcline84"> 84</a></span><span class="line">    <span class="keyword">end</span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,85" id="srcline85"> 85</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,86" id="srcline86"> 86</a></span><span class="line">    <span class="comment">% Compute the RHS for step of method</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,87" id="srcline87"> 87</a></span><span class="line">    <span class="mxinfo" id="T9:U115"><span class="mxinfo" id="T9:U116"><span class="var type1" id="S26T11U578">f</span>(:,<span class="mxinfo" id="T1:U118">1</span>)</span> = <span class="mxinfo" id="T9:U119"><span class="fcn" id="F115N2:B582">fy_</span>(<span class="var type1" id="S18T1U583">t</span> , <span class="var type1" id="S25T9U584">x</span> ,  <span class="var type1" id="S7T4U585">auxdata</span>)</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,88" id="srcline88"> 88</a></span><span class="line">    <span class="keyword">for</span> <span class="var type1" id="S29T1U588">j</span> = <span class="mxinfo" id="T1:U124">1</span>: 12,</span></span>
<span class="srcline"><span class="lineno"><a href="1,89" id="srcline89"> 89</a></span><span class="line">        <span class="mxinfo" id="T9:U125"><span class="mxinfo" id="T9:U126"><span class="var type1" id="S26T11U595">f</span>(:,<span class="mxinfo" id="T1:U128"><span class="var type1" id="S29T1U598">j</span>+<span class="mxinfo" id="T1:U130">1</span></span>)</span> = <span class="mxinfo" id="T9:U131"><span class="fcn" id="F115N2:B601">fy_</span>(<span class="var type1" id="S18T1U603">t</span>+<span class="mxinfo" id="T1:U133"><span class="mxinfo" id="T1:U134"><span class="var type1" id="S8T26U606">c_i</span>(<span class="var type1" id="S29T1U607">j</span>)</span>*<span class="var type1" id="S14T1U608">h</span></span>, <span class="mxinfo" id="T9:U138"><span class="var type1" id="S25T9U610">x</span>+<span class="mxinfo" id="T9:U140"><span class="mxinfo" id="T11:U141"><span class="var type1" id="S14T1U613">h</span>*<span class="var type1" id="S26T11U614">f</span></span>*<span class="mxinfo" id="T28:U144"><span class="var type1" id="S9T27U616">a_i_j</span>(:,<span class="var type1" id="S29T1U618">j</span>)</span></span></span>, <span class="var type1" id="S7T4U619">auxdata</span>)</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,90" id="srcline90"> 90</a></span><span class="line">    <span class="keyword">end</span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,91" id="srcline91"> 91</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,92" id="srcline92"> 92</a></span><span class="line">    <span class="comment">%%%%%%%%%% Forced exit condition %%%%%%%%%%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,93" id="srcline93"> 93</a></span><span class="line">    <span class="mxinfo" id="T1:U148"><span class="var type1" id="S24T1U622">fk</span> = <span class="mxinfo" id="T1:U150"><span class="var type1" id="S24T1U624">fk</span> + <span class="mxinfo" id="T1:U152">13</span></span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,94" id="srcline94"> 94</a></span><span class="line">    <span class="keyword">if</span> <span class="mxinfo" id="T2:U153"><span class="var type1" id="S24T1U629">fk</span> &gt; <span class="mxinfo" id="T1:U155">5e4</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,95" id="srcline95"> 95</a></span><span class="line">        <span class="mxinfo" id="T1:U156"><span class="var type1" id="S4T1U633">ErrorFlag</span> = <span class="mxinfo" id="T1:U158">1</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,96" id="srcline96"> 96</a></span><span class="line">        <span class="keyword">return</span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,97" id="srcline97"> 97</a></span><span class="line">    <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,98" id="srcline98"> 98</a></span><span class="line">    <span class="comment">%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,99" id="srcline99"> 99</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,100" id="srcline100">100</a></span><span class="line">    <span class="comment">% Two solution</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,101" id="srcline101">101</a></span><span class="line">    <span class="mxinfo" id="T9:U159"><span class="var type1" id="S30T9U638">sol2</span> = <span class="mxinfo" id="T9:U161"><span class="var type1" id="S25T9U640">x</span> + <span class="mxinfo" id="T9:U163"><span class="mxinfo" id="T11:U164"><span class="var type1" id="S14T1U643">h</span>*<span class="var type1" id="S26T11U644">f</span></span>*<span class="mxinfo" id="T28:U167"><span class="var type1" id="S10T28U645">b_8</span></span></span></span></span> ; <span class="comment">% 8-th order solution</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,102" id="srcline102">102</a></span><span class="line">    <span class="mxinfo" id="T9:U169"><span class="var type1" id="S31T9U648">sol1</span> = <span class="mxinfo" id="T9:U171"><span class="var type1" id="S25T9U650">x</span> + <span class="mxinfo" id="T9:U173"><span class="mxinfo" id="T11:U174"><span class="var type1" id="S14T1U653">h</span>*<span class="var type1" id="S26T11U654">f</span></span>*<span class="mxinfo" id="T28:U177"><span class="var type1" id="S11T28U655">b_7</span></span></span></span></span> ; <span class="comment">% 7-th order solution</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,103" id="srcline103">103</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,104" id="srcline104">104</a></span><span class="line">    <span class="comment">% Truncation error</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,105" id="srcline105">105</a></span><span class="line">    <span class="mxinfo" id="T1:U179"><span class="var type1" id="S32T1U658">error_1</span> = <span class="mxinfo" id="T1:U181">norm(<span class="mxinfo" id="T9:U182"><span class="var type1" id="S31T9U662">sol1</span>-<span class="var type1" id="S30T9U663">sol2</span></span>)</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,106" id="srcline106">106</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,107" id="srcline107">107</a></span><span class="line">    <span class="comment">% Estimate the error and the acceptable error</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,108" id="srcline108">108</a></span><span class="line">    <span class="mxinfo" id="T1:U185"><span class="var type1" id="S34T1U666">Error_step</span> = <span class="mxinfo" id="T1:U187">norm(<span class="var type1" id="S32T1U669">error_1</span>,<span class="string">'inf'</span>)</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,109" id="srcline109">109</a></span><span class="line">    <span class="mxinfo" id="T1:U189"><span class="var type1" id="S35T1U673">tau</span> = <span class="mxinfo" id="T1:U191"><span class="var type1" id="S15T1U675">tol</span>*<span class="mxinfo" id="T1:U193">max(<span class="mxinfo" id="T1:U194">norm(<span class="var type1" id="S25T9U680">x</span>,<span class="string">'inf'</span>)</span>,<span class="mxinfo" id="T1:U196">1.0</span>)</span></span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,110" id="srcline110">110</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,111" id="srcline111">111</a></span><span class="line">    <span class="comment">% Update the solution only if the error is acceptable</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,112" id="srcline112">112</a></span><span class="line">    <span class="keyword">if</span> <span class="mxinfo" id="T2:U197"><span class="var type1" id="S34T1U686">Error_step</span> &lt;= <span class="var type1" id="S35T1U687">tau</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,113" id="srcline113">113</a></span><span class="line">        <span class="mxinfo" id="T1:U200"><span class="var type1" id="S18T1U690">t</span> = <span class="mxinfo" id="T1:U202"><span class="var type1" id="S18T1U692">t</span> + <span class="var type1" id="S14T1U693">h</span></span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,114" id="srcline114">114</a></span><span class="line">        <span class="mxinfo" id="T9:U205"><span class="var type1" id="S25T9U696">x</span> = <span class="var type1" id="S30T9U697">sol2</span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,115" id="srcline115">115</a></span><span class="line">        </span></span>
<span class="srcline"><span class="lineno"><a href="1,116" id="srcline116">116</a></span><span class="line">        <span class="mxinfo" id="T1:U208"><span class="var type1" id="S2T1U700">tout</span> = <span class="var type1" id="S18T1U701">t</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,117" id="srcline117">117</a></span><span class="line">        <span class="mxinfo" id="T35:U211"><span class="var type1" id="S3T35U704">xout</span> = <span class="mxinfo" id="T35:U213"><span class="var type1" id="S25T9U706">x</span>.'</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,118" id="srcline118">118</a></span><span class="line">        </span></span>
<span class="srcline"><span class="lineno"><a href="1,119" id="srcline119">119</a></span><span class="line">        <span class="mxinfo" id="T1:U215"><span class="var type1" id="S23T1U709">reject</span> = <span class="mxinfo" id="T1:U217"><span class="var type1" id="S23T1U711">reject</span> - <span class="mxinfo" id="T1:U219">1</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,120" id="srcline120">120</a></span><span class="line">    <span class="keyword">else</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,121" id="srcline121">121</a></span><span class="line">        <span class="mxinfo" id="T1:U220"><span class="var type1" id="S22T1U716">n_reject</span> = <span class="mxinfo" id="T1:U222"><span class="var type1" id="S22T1U718">n_reject</span> + <span class="mxinfo" id="T1:U224">1</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,122" id="srcline122">122</a></span><span class="line">        <span class="mxinfo" id="T1:U225"><span class="var type1" id="S23T1U722">reject</span> = <span class="mxinfo" id="T1:U227">1</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,123" id="srcline123">123</a></span><span class="line">    <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,124" id="srcline124">124</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,125" id="srcline125">125</a></span><span class="line">    <span class="comment">% Step control</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,126" id="srcline126">126</a></span><span class="line">    <span class="keyword">if</span> <span class="mxinfo" id="T2:U228"><span class="var type1" id="S34T1U727">Error_step</span> == <span class="mxinfo" id="T1:U230">0.0</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,127" id="srcline127">127</a></span><span class="line">        <span class="mxinfo" id="T1:U231"><span class="var type1" id="S34T1U731">Error_step</span> = <span class="mxinfo" id="T1:U233"><span class="mxinfo" id="T1:U234">eps</span>*<span class="mxinfo" id="T1:U235">10.0</span></span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,128" id="srcline128">128</a></span><span class="line">    <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,129" id="srcline129">129</a></span><span class="line">    <span class="mxinfo" id="T1:U236"><span class="var type1" id="S14T1U738">h</span> = <span class="mxinfo" id="T1:U238">min(<span class="var type1" id="S13T1U741">hmax</span>, <span class="mxinfo" id="T1:U240"><span class="mxinfo" id="T1:U241"><span class="mxinfo" id="T1:U242">0.9</span>*<span class="var type1" id="S14T1U745">h</span></span>*<span class="mxinfo" id="T1:U244">(<span class="mxinfo" id="T1:U245"><span class="var type1" id="S35T1U749">tau</span>/<span class="var type1" id="S34T1U750">Error_step</span></span>)^<span class="var type1" id="S12T1U751">pow</span></span></span>)</span></span>;</span></span>
<span class="srcline"><span class="lineno"><a href="1,130" id="srcline130">130</a></span><span class="line">    <span class="keyword">if</span> (<span class="mxinfo" id="T2:U249"><span class="mxinfo" id="T1:U250">abs(<span class="var type1" id="S14T1U758">h</span>)</span> &lt;= <span class="mxinfo" id="T1:U252">eps</span></span>)</span></span>
<span class="srcline"><span class="lineno"><a href="1,131" id="srcline131">131</a></span><span class="line">        <span class="keyword">if</span> <span class="mxinfo" id="T2:U253"><span class="var type1" id="S23T1U764">reject</span> == <span class="mxinfo" id="T1:U255">0</span></span></span></span>
<span class="srcline"><span class="lineno"><a href="1,132" id="srcline132">132</a></span><span class="line">            <span class="comment">% warning</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,133" id="srcline133">133</a></span><span class="line">            <span class="comment">%             disp('Warning!!! ode87. Step is very small !!!') ;</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,134" id="srcline134">134</a></span><span class="line">            <span class="mxinfo" id="T1:U256"><span class="var type1" id="S14T1U768">h</span> = <span class="mxinfo" id="T1:U258"><span class="mxinfo" id="T1:U259">eps</span> * <span class="mxinfo" id="T1:U260">100</span></span></span> ;</span></span>
<span class="srcline"><span class="lineno"><a href="1,135" id="srcline135">135</a></span><span class="line">            <span class="keyword">return</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,136" id="srcline136">136</a></span><span class="line">        <span class="keyword">else</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,137" id="srcline137">137</a></span><span class="line">            <span class="comment">% error</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,138" id="srcline138">138</a></span><span class="line">            <span class="comment">%             disp('Error in ode87. Step is too small.') ;</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,139" id="srcline139">139</a></span><span class="line">            <span class="keyword">return</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,140" id="srcline140">140</a></span><span class="line">        <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,141" id="srcline141">141</a></span><span class="line">    <span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,142" id="srcline142">142</a></span><span class="line">    </span></span>
<span class="srcline"><span class="lineno"><a href="1,143" id="srcline143">143</a></span><span class="line"><span class="keyword">end</span></span></span>
<span class="srcline"><span class="lineno"><a href="1,144" id="srcline144">144</a></span><span class="line"></span></span>
<span class="srcline"><span class="lineno"><a href="1,145" id="srcline145">145</a></span><span class="line"><span class="keyword">end</span></span></span>
</pre>
